
#include "libcmat.h"

CMatrix::CMatrix(const char* refFileName, const char* dbsnpFileName) {
	_refFile.open(refFileName);
	if(dbsnpFileName) {
		_dbsnpFile.open(dbsnpFileName);
	}
	cout << "refFile: " << refFileName << endl;
	cout << "dbsnpFile: " << dbsnpFileName << endl;
	_genome = new Genome(_refFile, _dbsnpFile);
	memset(_countMatrix, 0, sizeof(ubit64_t)*cmatrixSize);
}


CMatrix::~CMatrix() {
	_refFile.close();
	_dbsnpFile.close();
	delete _genome;
}

void CMatrix::update(unsigned int readLen, char* read, char* qual, 
			unsigned int position, int hit, char strand, 
			const char* chrName) {
	map<Chr_name, Chr_info*>::iterator current_chr;
	current_chr = _genome->chromosomes.find(chrName);//_genome->chromosomes.end();
	ubit64_t ref(0);

	if (hit == 1) {
		for(unsigned int coord = 0; coord != readLen; coord++) {
			if (read[coord] == 'N' || read[coord] == 'n') {
				;
			}
			else {
				ref = current_chr->second->get_bin_base(position + coord);
				if ( (ref&12) !=0 ) {
					// This is an N on reference or a dbSNP
					;
				}
				else {
					if(strand=='+') {
						// forward strand
						_countMatrix[(((ubit64_t)qual[coord])<<12) | (coord<<4) | ((ref&0x3)<<2) | ((read[coord]>>1)&3)] += 1;
					}
					else {
						// reverse strand
						_countMatrix[(((ubit64_t)qual[coord])<<12) | ((readLen-1-coord)<<4) | ((ref&0x3)<<2) | ((read[coord]>>1)&3)] += 1;
					}
				}
			}
		}
	}
}

void CMatrix::write(const char* fileName) {
	FILE* file = NULL;
	safe_fopen(&file, fileName, "wb");
	safe_fwrite(_countMatrix, sizeof(ubit64_t), cmatrixSize, file);
	safe_fclose(file);
}

